
'''
Created on Jul 2, 2011

@author: oabalbin

It calculates the coverage over the whole or particular positions
of the genome.
Uses samtools pileup and varscan readcounts for doing that
'''

import subprocess

def calculate_read_coverage(bam_file, ref_genome, path_to_sam, path_to_varscan,
                            output_file, target_regions=None):
    '''
    Calculates a pileup using or not target regions
    Uses varscan to calculate the read counts a particular site.
    '''
    min_coverage=str(1)
    min_base_quality=str(20)
    samtools_command=path_to_sam+'samtools'
    varscan_command = path_to_varscan+'varscan'
    if target_regions is not None:
        args=[samtools_command,'pileup','-s','-f',ref_genome] + [bam_file] +['-l', target_regions]+\
            ['|', 'java','-jar',varscan_command, 'readcounts', '--min-coverage', min_coverage, '--min-base-qual',
             min_base_quality, '--output-file',output_file]
    else:
        args=[samtools_command,'pileup','-s','-f',ref_genome] + [bam_file] + \
            ['|', 'java','-jar',varscan_command, 'readcounts', '--min-coverage', min_coverage, '--min-base-qual',
             min_base_quality, '--output-file',output_file]
        
    args= [a.replace(',',';') for a in args]
    args = ",".join(args).replace(',',' ').replace(';',',')
    
    return args



if __name__ == '__main__':
    
    path_to_sam = '/exds/users/oabalbin/sw/samtools-0.1.10/' 
    ref_genome = '/exds/projects/alignment_indexes/gatk/hg19/hg19.fa'
    bam_file = '/exds/users/oabalbin/projects/exomes/HCC1599/sample/HCC1599_n/HCC1599_n_merged.quick.smdup.bam' 
    target_regions= '/exds/users/oabalbin/projects/exomes/genelists/informative_genes_cosmic_mutations.txt.ActionableGeneCosmic.pos' 
    path_to_varscan = '/exds/users/oabalbin/sw/SNVcallers/'
    output_file = '/exds/users/oabalbin/projects/exomes/HCC1599/sample/HCC1599_n/HCC1599_n_merged.quick.actionable.pileup'
    
    args = calculate_read_coverage(bam_file, ref_genome, path_to_sam, path_to_varscan,
                            output_file, target_regions)
    
    print args
    
    ret = subprocess.call(args)



